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We propose a finite-difference algoritlrm for solving tfie time-dependent Ginzburg- 
Landau (TDGL) equation coupled to the appropriate Maxwell equation. The time 
derivatives are discretized using a second order semi-implicit scheme which, for 
intermediate values of the Ginzburg-Landau parameter k, allows time-steps two 
orders of magnitude larger than commonly used in explicit schemes. We demonstrate 
the use of the method by solving a fully three-dimensional problem of a current- 
carrying wire with longitudinal and transverse magnetic fields. © ??? Academic Press 
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1. GINZBURG-LANDAU MODEL 

In the Ginzburg-Landau model, a superconductor is characterised by a complex order 
parameter ijj. The local density of superconducting electrons is represented by ji/ip. The 
theory postulates that close to the critical temperature, the free energy can be expanded in 
a series of the form 
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where a and b are phenomenological parameters that depend on external parameters such 
as temperature, A denotes the vector potential, H an external magnetic field, and eg and 
rua are the effective charge and the effective mass of the Cooper pairs. Below the transition 
temperature T^, a becomes negative, whereas 6 > for all T. 

1.1. The time-dependent Ginzburg-Landau (TDGL) equations 

The equations of motion for the order parameter and the vector potential are the Euler- 
Lagrange equations of the free energy functional, 

+ = ^(v-i|A)%+|a|V'-6|V'lV (2) 
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j„ = C7(-V$ - dtA) , (5) 

where £) is a phenomenological diffusion constant, and $ is the electric potential, included 
to retain the gauge invariance of the equations. Equation (3) is the Maxwell equation for the 
magnetic field, where the displacement current cq-B has been omitted as it only becomes 
significant for velocities close to the speed of fight. The total current is given by the sum 
of the supercurrent, j^, and the normal current, j„, which obeys Ohm's law. 

1.2. Dimensionless units 

We scale length in multiples of the coherence length, ^ = h/ ^/2m\a\, time in r = 
the wavefunction in (/;() = ^\a\/b, the vector potential in Aq = \/2kH where 
He = /xo|op/fo, the electric potential in $o = (C/''")^o and resistivity in units of the 
normal resistivity (To = I/k^-D/xo- The so-caUed Ginzburg-Landau parameter is given by 
= 2171% / e^T^ jjLQ. The characteristic length scale for variations of the magnetic field is 
A = and V x measures the magnetic field in units of ^/2kHc = -ffc2- In scaled units 
equations (2) and (3) become 

+ = (V-iA)^^ + V- IV-IV (6) 

K^v X V X A = (VS - + (-V$ - dtA) + k^V x if , (7) 
^ , ' ^ , ' ' — ^ ' 

3. 3n Je»* 

where S denotes the phase of tj^j. The last term in equation (7) can be understood as 
an external current jg^.^ with Vj^xt = 0- the following, this term will be omitted. 
However, it can be easily included in the algorithm, for instance, to model magnetic 
impurities. In dimensionless units, the dynamics of the superconductor depends on the 
dimensionless Ginzburg-Landau parameter k only. For values k < l/\/2 one finds a 
behaviour characteristic of a type-1 superconductor whereas for k > 1/^/2 a type-II 
superconductor is modelled. 

1.3. Gauge transformation 

The dynamics of the measurable quantities E, B, and j are invariant under the 
transformation 



(8) 



where A is an arbitrary scalar field. We choose the zero potential gauge, A(r,t) = 
/ dt^{r, t), in other words, $(r) = at all times. For this choice, equations (6) and (7) 
become 

dti^ = (V - iA)^ 7/; + - IV'PV^ (9) 
dtA = (VS'-A)IV'P - K^V x V x A . (10) 
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In the following section, we suggest a fast and reliable numerical method to find an 
approximate solution to these equations. 
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2. NUMEMCAL METHODS 

The most popular approach to the solution of the TDGL equations, (9) and (10), is a 
gauge-invariant discretization that is second order accurate in space and first order in time 
[1, 2, 3, 4, 5]. In addition, a number of other finite difference [6, 7] and finite element 
methods [8, 9] have been developed. For large values of k, the magnetic field is nearly 
homogeneous and equation (10) can be dropped. This case is often referred to as the London 
limit. The remaining equation has been solved by a semi-implicit Fourier spectral method 
which is second order accurate in time [10]. An equation very similar to equation (9), 
the Gross-Pitaevskii equation, is used to model vortex dynamics in dilute Bose-Einstein 
condensates [11]. Here, we modify the very robust and accurate semi-implicit Crank- 
Nicholson algorithm used in [1 1] to include the equation for the vector potential. 



2.1. The U -ip method 

The widely used U — tjj method is described in detail by Gropp et. al. [1]. As this 
method forms the basis of our algorithm we briefly review the main points here. Complex 
link variables U^, and are introduced to preserve the gauge invariant properties of 
the discretized equations. 



U''{x,y,z) = exp l~i {x' ,y, z)dx' 

\ Jxo y 

Uy{x, y, z) = exp (^-i j"^ Ay{x, y' , z)dy'^ (11) 
U''{x,y,z) =exp(^-i A''{x,y, z')dz'^ , 



where {xo,yo,zo) is an arbitrary reference point. The TDGL equations can then be 
expressed as functions of t/j and these link variables. Both the order parameter and the 
link variables are discretized on a three dimensional grid with grid spacing hx, hy, and 
hz, respectively. The mesh points for the link variables are half way between the mesh 
points for the order parameter (see Fig. 1). AH spatial derivatives are approximated by 
finite differences to second order accuracy. Denoting the complex conjugate of U by U, 
the finite difference representations of the TDGL equations read 



dti>i,j,k 



+ 



+ 



hi 

^,j-i,fcV'i,3-i,fc - 2V'»,j,fc + Ulj,k'^i,j+i,k 

fly 

ulj,k-i'^i,j,k-i - "^i^ij.k + Ui,j,k'^i,j,k+i 



hi 



-h (1 - \tpij,k\'^)->pi,j,k 



dtUrj,k = -ilni(^i:,,fc)t/i:,-,fc, 



(12) 
(13) 
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where 

rx _ 2 '.i+l.fe '.i.fc ^ i,3,k'^ i,j-l,k^i,j-l,k^i+l,j-l,k 

•^i,j,k - 1^ L2 

—Z —X —Z —X 

2 ^ i+l,3,k-lU i,j,k-lUi,j,k-l^i,j,k ~ ^ i+l,j,kU i,3,k^i,j,k^i,j,k+l 
+ U^,j,k^i,j,k^i+h3,k ■ 

Analogous expressions for dtUV - ^. and dtU^ ^ ^. can be obtained by permutating the coor- 
dinates and indices as follows, 

{x, y, z; i, j, k) {y, z, x; j, k, i) {z, x, y; k, i, j) {x, y, z; i, j, k) . (14) 

The time evolution is approximated by a simple Euler step, 

V-^-feCt + At) = ^Pi^j^kit) + At dt^^^j^kit) + OiAf) (15) 
f/;^^. ,(t + At) = C/-^. + At dtU^^^^,{t) + 0{At^) . (16) 

To keep U^jj. uni-modular, equation (16) is often modified to 

Urj,k{t + Ai) = Urj,kit) exp i-iAt lmj^l3,k) + OiAt^) • (17) 

The Euler method is only first order accurate in time, i.e., the truncation error made due to 
the finite difference approximation of the time derivative is proportional to Af^. However, 
the main problem is that the code becomes unstable if long time steps are used. The cause 
of this instability is the diffusion-like character of the dynamics described by the equations 
(12) and (13). Equation (12) can innmediately be written as a diffusion equation with an 
additional non-linear term 

dt'4'i,j,k = Lx'4'i,j,k + Lyipij^k + Lzi}i,j,k + f , (18) 

where / stands for (1 — | V'i,j,fe |^)V'i,j,fc -^x> Ly, and denote the weighted Laplacian 
operators, 

J- , _ ai-li>i-l,j,k - '^i'ij.k + ai+l'>Pi+l,j,k ,1„^ 

i^xVi,j,k = , (ly) 

with |ai_i| = \ai+i\ = 1 in our case. The diffusion constant is 1 in dimensionless units. 
Equation (13) is also dominated by diffusive terms as will become evident in the next 
section. The diffusion constant for the vector potential is k^. This can be seen by taking 
the curl of equation (10), B = k^V'^B -|- V x j^. The one-step forward Euler method is 
only stable as long as the time step is shorter than the diffusion time across a cell of width 
/i [13]. For example, using a grid spacing of /i = 0.5^ and k = 4, the theoretical limit for 
the time step is 

^*<iJ = S"«-«°^^- ^^^^ 

In practice, a time step of At = 0.0025 is used to ensure stability [1]. In contrast, a semi- 
implicit two-step algorithm is unconditionally stable for diffusive problems and enables 
much larger time-steps to be employed. 
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2.2. Semi-implicit algoritlmi 

We propose a spatial discretization of the equations very similar to the above U — ip 
method. The link variables are uni-modular, \U^ji^\ = 1, and can be written as the 
exponential of a phase, U^ jj. = exp{—i(|)^ J ,,). We use the real- valued variable instead 
of the complex- valued W. The fields tp and <p are represented on a three-dimensional 
grid. The mesh points of the phase factors are placed between the mesh points of the order 
parameter (see Fig. 1). For the field V'j.j.fe^ the grid point indices are i ~ l...Nx + 1, 
j = l...Ny + 1, and k = l.-.N^ + 1. For the indices in the x direction range 

i = l...Nx only, due to the relative displacement of the grids. Similarly, j = l...Ny for 
Cfc^dfc = 1...7V,for<^f_.^,. 



\|/ d) \|/ d) \|/ — 

i-i,j+i 



^' ^' 



\|/ d) \|/ d) \|/ — 



B^ B^ 



\|/ d) \|/ d) \|/ — 



FIG. 1. The evaluation points for the fields f/' and </> in the a; — y plane. A finite difference approximation 
for the magnetic field is given in (36). 

We now discretize the spatial derivatives in equations (9) and (10) using the modified 
link variables (fr^ , (ffl, and 0^. For equation (9), we can reuse the expansion (12) except 
that U^j^f, is replaced by exp(— i(^^_^^^), etc. 



+ 



+ 



exp(i<^?_ij^fc)V'»-i,j,fc - 2V>j,j- fc -I- exp{-i(f)lj ^.)tpi+ij^k 

hi 

exp(i?!'fj_i,fe)V'i,j-i,fc - '2ipi,j,k + exp(-i^*'^^._^)V'i,j+i,fe 

"'y 



(21) 



+ (1 - \fpi,j,k\'^)'ipi,j,k ■ 



With help of the relation — VxVxA = V^A — V{VA), the second order accurate finite 
difference representation of (10) is 



6 



T. WINIECKI AND C. S. ADAMS 



'V 



(22) 



The expressions for dtcj)^ j ^ and dt4>ij i. given by cyclic permutation (14). 
Note, that the discretized equations are still invariant under the gauge transformation 



4>lj,k €,j^k + (Aij+i,fe - A,j- fe) 
•^fj.fc ^ (l'lj,k + (A»j,fc+i - A,,j- fc) 



(23) 



Retaining the gauge invariance at the discrete level is often equivalent to preserving certain 
conservation laws and physical principles. It is crucial that the numerical approximation 
does not depend on the particular choice of gauge. If, for example, one studies the motion 
of a vortex lattice due to an applied electric field E^, the measurable quantities B, 
and j oscillate in time [16]. The system is driven through a series of equivalent solutions 
and the dynamics is roughly described by A = E^xt. This means that the phase gradients 
in the order parameter build up in time and the phase difference between two neighbouring 
grid points eventually exceeds 27r. This is normally a problem as the finite difference 
approximation becomes invalid. However, using the link variables U or (p these phase 
gradients are exactly cancelled by the change in the vector potential. 

We now want to introduce a new scheme to update the wavefunction ?/'*^"^and the link 
variables cj)^"^ from the nth to the (n + l)th step. The idea is to treat the diffusive terms 
semi-impUcitly whereas all other terms are still treated explicitly. In this way we reduce 
the stability constraints associated with the simple Euler method but avoid the expensive 
solution of non-linear equations. The technique is known as the method of fractional steps 
[12]. A second-order accuracy in the time-step can be achieved by a simple 3 step iteration. 

As mentioned above, the diffusive character of equation (10) becomes apparent in the 
new discretization and both equation (21) and equation (22) can be written as a initial value 
problem of the form 



where u stands for the fields tp or (p^, (p^ or (p^, respectively, D is the diffusion constant 
with D ~ in (22) and / indicates all the other terms, (1 — |'i/'j,j,fe|^)V'i,j,fe in (21) and 
the last three fines in equation (22). Note, that i^; = in (22). 

The second derivatives are approximated in the usual way by an expression involv- 
ing three neighbouring grid points. For any pair (j, k) the action of on the vector 
{ui,j,k}i=2...Na,, can be represented by a tri-diagonal matrix S^, 



dtUij^k = D {LxUij^k + LyUij^k + LzUij^k) + f 



(24) 




'X 




(25) 
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As emphasized before, the instabilities of the Euler method have their origin in the expUcit 
treatment of the diffusive terms. We now discretize the time derivative in equation (24) in 
the following way, 



At 2 hi 

ID51 



2 m 



This discretization is semi-implicit as the right hand side of the equation depends on the 
fields at the old and the new time level. This mixing leads to an improved accuracy and 
prevents the algorithm from developing instabihties. After rearranging the equation we get 



2hl ^ 2hl y 2hl ^ ' 



(27) 



We now employ an approximate factorisation [12], 



1 _ :^a2 _ :^x2 _ 2^x2 1 ^ 
2hl ^ 2hl y 2hl ^ ' 

, DAt / DAt DAt 



2kVVV 2.?-^^ ^-W'^^)+^^^*^^• 



(28) 



The multidimensional operator is spht into three operators that involve difference approxi- 
mations in only one dimension. With the abbreviations 

DAt / DAt „\ 

and considering — u^")) = 0{At) equation (27) becomes 

= B^ByB.u^^^ + ^ + /(")) + 0{At^) . (30) 

The tri-diagonal matrices A and B are actually time dependent because the differential 
operators L in equation (21) depend on the link variables. In the above equation, A is a 
function of (/)("+^) whereas B depends on 0^"' . Consequently, the equations are solved in 
the following step-wise manner 

A.u("+V3) = B^ByB^u^-) + ^ (/("+!) + /(")) 

A^„("+2/3) = uin+l/S) (31) 

where the 'fractional' time levels indicate intermediate results. The explicit term as 
well as the matrix elements of A may depend on the values of the link variables at the new 
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time level and are unknown initially. We assume that = w^") to start with. After 

the first iteration of equation (31) for all variables tjj and (f), the updated values at the new 
time level are used in the matrix elements of A for the second iteration, and so on. The 
product BxByBzu'^"-^ is a function of known values at the previous time level and can be 
stored in an auxiUary variable for subsequent iterations. As the matrices A are tri-diagonal, 
fast inversion routines can be apphed [13]. 

The entire algorithm relies on the convergence of this iteration technique. To test whether 
the procedure converges, we calculate a total update, S{m), of the fields after m iterations 
by comparing all values at the time level (n + 1) to all values at the previous time step, (n), 



Sim) = ^(( 

i,j,k \ ^ 



, (32) 



where the three dots indicate the corresponding terms for the fields and Fig. 2 
shows a typical evolution of the update for a time step of At = 0.5. After as few as five 
iterations the approximated increment is very close to the exact value. For smaller time 
steps, the procedure converges faster. We find an optimum trade-off between accuracy and 
performance for three iterations. We further check, if the correction between two successive 
iterations. 



T{m) 




(n+l,m) 



, (n+l,m-l) 



;x(n+l,m— 1) 



(33) 



converges to zero. Fig. (2, inset) confirms an exponential convergence. 




5 10 15 20 

iteration 



FIG. 2. The modification to the solution after one time step (total update) versus the number of iterations. For 

smaller time-steps, the solution converges faster We find an optimum of speed and accuracy for a combination 
of three iterations and a time-step of At = 0.5. The correction between iterations, T(m) is shown inset. 
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The accuracy of the method is assessed by comparing the solution to simulations using 
the Euler method with a much smaller time step, At = 0.0025. Up to a time step of 
At = 0.5, no significant deviations could be observed. The program runs about 40 times 
faster than the Euler method for these parameters. For the calculations below we use a finer 
grid (h = 0.4) and a sUghtly larger Ginzburg-Landau parameter (k = 5). The speed-up for 
these values is about 100. 

Our implicit method is less memory intensive than the standard U — tp method because 
it uses real- valued link variables rather than complex- valued ones that must be represented 
by two real numbers. For a grid of points, the Euler method uses an equivalent of 22N^ 
real-valued variables (the complex wavefunction and the three complex link variables at 
two time levels plus three complex fields W, see [1]) whereas the imphcit method requires 
the storage of a total of 19 variables (the complex wavefunction and the three real Unk 
variables at two time levels, the products B^ByBz u(") in (31) plus four auxiUary fields). 



The correct implementation of the boundary conditions requires great care because of the 

relative displacement of the grids. The matrices and only act on the interior points, 
i = 2...Nx, of the vectors V'i, (j^i, and^f for all j = l...Ny + l,k = l...N;, + l. Note, that 
Ax = Bx = 1 in the case u = c/)^. The end points i = 1 and i = N^ + laie computed for 
book-keeping purposes. Similarly, the operators By and B^ do not automatically include 
information on the end points at j = 1, Ny + 1 and k = 1 , TV^ + 1 , respectively. In addition, 
there are different boundary conditions, namely periodic, Dirichlet and Neumann boundary 
conditions, that require an adaption of the matrix elements on the first and last row [13]. 
Another complication is that the physical boundary conditions that apply for the vectors 
and in equation (31) do not necessarily apply for the intermediate results 

^(n+i/3) gjj^ ^(n+2/3) ^^hcn, for example, solving the second equation of the system 



a correct treatment of the boundary condition for = ^^17/"+^^ must be imple- 

mented into Ay. The matrices A,., Ay, and A^ commute as they act in different directions. 
It is therefore advisable to solve (31) starting in the direction with the simplest boundary 
condition. For a periodic boundary condition, for example, the relations u\ = Un and 
Ujv+i = 1*2 hold at all time levels, including 'fractional' ones. 

The boundary conditions depend on the geometry of the problem. We choose a system 
with a periodic boundary conditions in the i;-direction. At the interfaces in the x- and 
the y-direction, boundary conditions for the magnetic field and the order parameter are 
applied. For the order parameter, V'i,i,fe. conditions are needed for all values at the faces 
of the three-dimensional box. The grid representation of the periodic boundary condition 
reads 



2.3. Boundary conditions 



(31), 



A^t,("+V3) = ^(-+1/3) 



(31') 




(34) 
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In the y and z direction we set the supercurrent across the boundary to zero [1], i. e., 

'>Pi,Ny + l,k 

Expressions for the end points of the link variables can be found by incorporating infor- 
mation of the magnetic field at the boundaries of the box. The three components of the 
magnetic field are given by the following second order finite-difference approximations. 

TJX 

By.. 

From these expressions, appropriate boundary conditions can be obtained. For example, 
the field (j^f^^i^ is unknown at j — 1, and we use the last equation to relate the values of 
(f>iik^^ known values 

<l^li,k = -BU,uh.K + 0^,2,. + <t^li,k - <t^l+i^,k ■ (37) 

Equation (21), (22) and (31) combined with the boundary conditions (34), (35), and (37) 
provide all the information needed to solve a three-dimensional problem. 

3. EXAMPLE 1: WIRE WITH LONGITUDINAL FIELD 

As an example, we model an infinite cylindrical wire in three dimensions with an external 
magnetic field, Hz, appUed along its axis. Such a configuration has been studied in [7, 14]. 
Experiments have shown that the magnetic field can increase the critical current down the 
wire [15]. 

In a type-11 superconductor, a sufficiently strong magnetic field, H^, will enter the 
cylinder and become trapped in vortex tubes aligned parallel to the axis of the cylinder. 
A current I along the wire induces an additional circular field H^{r) such that I = 
2'KpK^H^{p) at a distance p from the axis of the wire. The vortex tubes corresponding to 
the current- induced field are rings coaxial with the wire. Above a critical current, these 
vortex rings enter at the edge of the cylinder and shrink until they annihilate on the axis of 
the cylinder. This process repeats and leads to dissipation. There is no stable mixed state 
associated with a field unless the vortex rings are pinned by impurities in the material. 
Blackburn et. al. [7] have argued, that by entangling the vortex rings with vortex lines due to 
a strong longitudinal field, the rings can be prevented from shrinking and thereby increase 
the critical current in the wire. 

We model a cylindrical shape of the wire in a rectangular box by adding a potential term 
y V with F = 5 to equation (9) at all grid points outside a cyUndrical region with radius 
R = 12^. The density of the order parameter outside the cylinder decreases rapidly to 
zero. An array of longitudinal vortex tubes is created by imposing boundary conditions 
for an external magnetic parallel to the wire. A current is ramped up by slowly increasing 



^/'2j,fcexp(-i(/.f_^.fe) 
V'i,2,feexp(-i(?if_i;,) 



(35) 



hyK ^ 



'^i,j,k+l 



'^i+l,j,k ■ 



^i,j,k + 4'i,j+l,k) 
^i^j,k 



(36) 
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FIG. 3. Time-independent arrangement of spiral vortices in a cylinder of radius 12^, k = 5. Five vortex 
lines are entangled with two rings. The applied fields are Hz = 0.2 and Hip = 1.5/ p, respectively. The tubes 
show the density of the order at the level \ip\'^ = 0.3. The left slice shows the magnetic field component parallel 
to the wire, H^, dark regions indicate a high field. The black lines mark the boundaries of the wire. 

a circular field (Bx and By) around the box until vortices enter. The Bean-Livingston 
surface barrier was lowered by adding a weak sinusoidal potential at the surface of the 
cylinder. Fig. 3 shows a time-independent state that arises after two vortex rings have 
entered the cylinder and entangled with the vortex lines. The critical current is dominated 
by the surface barrier as expected for small samples. Consequently, we do not observe any 
improvement in the critical current due to the presence of a longitudinal field. However, 
the effect could become more significant for larger sample sizes or if the surface effects are 
suppressed [16]. 

4. EXAMPLE 2: WIRE WITH TRANSVERSE FIELD 

In superconducting magnets, the external magnetic field is typically aligned perpendic- 
ular to the wire (B^ for example). In the mixed state, an array of vortex lines fills up the 
superconductor (see Fig. (4)). Any current carried by the wire superimposes a circular 
field onto this appUed field. As a result, a gradient of the magnetic field develops that 
can be associated with a Lorentz force on the vortices. In most appUcations, different 
pinning mechanisms balance this Lorentz force and freeze the flux lattice up to a critical 
current density. For larger currents, the Lorentz force exceeds the pinning force and vor- 
tices start moving [16]. The motion of the flux lattice coincides with the breakdown of 
superconductivity. 

In this geometry the magnetic field at the boundary of the computational box strongly 
depends on the currents inside unless the box size is much larger than the radius R of 
the cylinder. In the Meissner state, for example, no flux lines penetrate the wire and 
supercurrents in the surface of the superconductor cancel the external field in the bulk of 
the wire. The field fines of these surface currents also extend outside the sample at length 
scales of order R. In the mixed state, the induced field is smaller and can be regarded 
as a small correction. To find a self-consistent solution, the fields induced by both the 
supercurrents and normal currents have to be added to the applied field and included in the 
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FIG. 4. Array of vortices in a wire exposed to a perpendicular magnetic field Bx = 0.4. The vortices tubes 
enter end exit the surface of the superconductor normally. The Ginzburg-Landau parameter is k = 5, the radius 
of the wire R = 8^. 

boundary conditions at each time step. We calculate the induced field iif ind('') using the 
Biot-Savart law, which in our units has the form, 



This calculation is computationally expensive. The integral is approximated by summing 
over all grid points for each boundary point requiring a total of O {N^ ) calculations, whereas 
a time step takes 0{N^) calculations for a box of grid points. With periodic boundary 
conditions, the integration must also be extended to regions outside the box. The effect of 
including iJind is to bend the vortex lines, especially near the top and the bottom of the 
sample as apparent in Fig 4. To model the motion of flux lines, a fully self-consistent time- 
dependent solution can be found by iterating the boundary conditions at each time-step. 
However, in practice the study of the motion of vortices above the critical current does not 
seem to be feasible. Possible ways around this problem are to increase the box size so that 
the induced currents can be neglected, to cut-off the integral in (38) at a certain distance 
from the boundary, | r — r' | < 7?., or to update H-md {f) only in larger intervals rather than 
every time step. The latter approach is especially suitable for finding time-independent 
solutions. 



In summary, we have demonstrated the use of a semi-implicit finite-difference scheme 
to solve the three-dimensional time-dependent Ginzburg-Landau equations. The method 
converges if a Crank-Nicholson scheme is applied to all Laplacian-type terms while all other 
terms are treated explicitly. We iterate each time step to achieve second order accuracy 
in time. For intermediate values of the Ginzburg-Landau parameter k, the method is 
stable and accurate for time-steps two orders of magnitude larger than used in the standard 
expUcit schemes. If the magnetic field at the surface of the computational box is partly 
due to currents inside, a self-consistent solution can be found by iteration. However, a full 
Biot-Savart calculation of the induced magnetic field remains computationally expensive. 
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